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SUMMARY 

An approach for assessing the delamination propagation and growth 
capabilities in commercial finite element codes was developed and 
demonstrated for the Virtual Crack Closure Technique (VCCT) 
implementations in ABAQUS®. The Double Cantilever Beam (DCB) specimen 
was chosen as an example. First, benchmark results to assess delamination 
propagation capabilities under static loading were created using models 
simulating specimens with different delamination lengths. For each 
delamination length modeled, the load and displacement at the load point were 
monitored. The mixed-mode strain energy release rate components were 
calculated along the delamination front across the width of the specimen. A 
failure index was calculated by correlating the results with the mixed-mode 
failure criterion of the graphite/epoxy material. The calculated critical loads 
and critical displacements for delamination onset for each delamination length 
modeled were used as a benchmark. The load/displacement relationship 
computed during automatic propagation should closely match the benchmark 
case. Second, starting from an initially straight front, the delamination was 
allowed to propagate based on the algorithms implemented in the commercial 
finite element software. The load-displacement relationship obtained from the 
propagation analysis results and the benchmark results were compared. Good 
agreements could be achieved by selecting the appropriate input parameters, 
which were determined in an iterative procedure. 
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A benchmark case to assess the fatigue growth prediction capabilities was also 
created. In a first step, the number of cycles to delamination onset was 
calculated from the material data for mode I fatigue delamination growth onset. 
Then, the number of cycles during stable delamination growth was obtained 
incrementally from the material data for mode I fatigue delamination 
propagation. For the benchmark case, where the results for delamination onset 
and growth were combined, the delamination length was calculated for an 
increasing total number of load cycles. A delamination growth prediction 
analysis that accounts for delamination fatigue onset as well as stable growth 
should therefore yield results that closely resemble the benchmark curve fit. 
Starting from an initially straight front, the delamination was allowed to grow 
under cyclic loading in a finite element model. The number of cycles to 
delamination onset and the number of cycles during stable delamination growth 
for each growth increment were obtained from the analysis. Input control 
parameters were varied to study the effect on the computed delamination 
growth. The benchmarks enabled the selection of the appropriate input 
parameters that yielded good agreements between the results obtained from the 
automated analysis and the benchmark results. Overall, the results are 
encouraging but further assessment for mixed-mode delamination is required. 

1: INTRODUCTION 

The virtual crack closure technique (VCCT) (Rybicki and Kanninen, 1977; 
Krueger, 2004) is widely used for computing energy release rates based on 
results from finite element analyses and to supply the mode separation required 
when using a mixed-mode fracture criterion to predict delamination onset or 
growth. As new methods for analyzing composite delamination are 
incorporated in finite element codes, the need for benchmarking becomes 
important since each code requires specific input parameters unique to its 
implementation. Once the parameters have been identified, they may then be 
used with confidence to model delamination growth for more complex 
configurations. 

The current paper summarizes the development of benchmark examples used 
to assess the automated delamination propagation and delamination growth 
capabilities in commercial finite element codes (Krueger, 2008, 2010; Orifici 
and Krueger, 2010). The creation of benchmark cases based on the Double 
Cantilever Beam (DCB) specimen are shown step-by-step. The creation step is 
independent of the software used. After creating the benchmark, the approach 
is demonstrated for the automated analysis in the commercial finite element 
code ABAQUS®. Starting from an initially straight front, the delamination is 
allowed to grow based on the algorithms implemented into the commercial 
finite element software. Input control parameters are varied to study the effect 
on the computed delamination propagation and growth. Results obtained from 
the automated analysis should closely match the results created manually to 
obtain the benchmark example. 
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2: SPECIMEN AND MODEL DESCRIPTION 

For the current numerical investigation, the Double Cantilever Beam (DCB), as 
shown in Figure 1, was chosen. The DCB specimen is used to determine the 
mode I interlaminar fracture toughness, G/o Typical two- and three- 
dimensional finite element models of the DCB specimen are shown in 
Figures 2 and 3. Along the length, all models were divided into different 
sections with different mesh refinement. The specimen was modeled with six 
elements through the specimen thic kn ess. The plane of delamination was 
modeled as a discrete discontinuity in the center of the specimen. To create the 
discrete discontinuity, each model was created from separate meshes for the 
upper and lower part of the specimens with identical nodal point coordinates in 
the plane of delamination. Two surfaces (top and bottom surface) were created 
on the meshes as shown in Figure 3. Additionally, a node set was created to 
identify the intact (bonded nodes) region. The initial delamination front is 
highlighted in red. The modelling is discussed in detail in related reports 
(Krueger, 2008, 2010). 

3: CREATING A BENCHMARK SOLUTION FOR STATIC LOADING 


First, models simulating specimens with seven different delamination lengths 
(30 mm < ao < 40mm) were analyzed. For each delamination length modeled, 
the reaction loads P at the location of the applied displacement were calculated 
and plotted versus the applied opening displacement <5/ 2 as shown in Figure 4. 
The mode I strain energy release rate was calculated along the delamination 
front across the width of the specimen. The critical load, P cr it, when the failure 
index in the center of the specimen (y/B= 0) reaches unity {Gj/G c = 1), can be 
calculated based on the relationship between load P and the energy release rate 


P 2 dCp 
2 dA 


( 1 ) 


In equation (1), Cp is the compliance of the specimen and dA is the increase in 
surface area corresponding to an incremental increase in load or displacement 
at fracture. The critical load P cr u and critical displacement d crit /2 were 
calculated for each delamination length modeled 



and the results were included in the load/displacement plots as shown in 
Figure 4 (solid red circles). The results indicate that, with increasing 
delamination length, less load is required to extend the delamination. This 
means that the DCB specimen exhibits unstable delamination propagation 
under load control. Therefore, prescribed opening displacements 6/ 2 were 
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applied in the analysis instead of nodal point loads P to avoid problems with 
numerical stability of the analysis. It was assumed that the critical 
load/displacement results can be used as a benchmark. For the delamination 
propagation, therefore, the load/displacement results obtained from the model 
of a DCB specimen with an initially straight delamination of a = 30 mm length 
should correspond to the critical load/displacement path (solid red line) in 
Figure 4. 

4: RESULTS FROM AUTOMATED PROPAGATION ANALYSES 

A set of example results are shown in Figure 5 where the computed resultant 
force (load P) at the tip of the DCB specimen is plotted versus the applied 
crack tip opening (d/2) for different input parameters. To overcome the 
convergence problems, the methods implemented in ABAQUS® were used 
individually to study the effects. First, global stabilization was added to the 
analysis. For a stabilization factor of 2x1 O' 5 , the stiffness changed to almost 
infinity once the critical load was reached causing the load to increase sharply 
(plotted in red). The load increased until a point was reached where the 
delamination propagation started and the load gradually decreased following a 
saw tooth curve with local rising and declining segments. The gradual load 
decrease followed the same trend as the benchmark curve (in grey) but is 
shifted toward higher loads. For a stabilization factor of 2xl0" 6 (in green), the 
same saw tooth pattern was observed but the average curve was in good 
agreement with the benchmark result. Second, contact stabilization was added 
to the analysis. For all combinations of stabilization factors and release 
tolerances, a saw tooth pattern was observed, where the peak values were in 
good agreement with the benchmark result (example plotted in blue). Third, 
viscous regularization was added to the analysis to overcome convergence 
problems. Convergence could not be achieved over a wide range of viscosity 
coefficients when a default release tolerance value of 0.2 was used. 
Subsequently, the release tolerance value was increased. For all combinations 
of the viscosity coefficient and release tolerance, where convergence was 
achieved, a saw tooth pattern was obtained, where the peak values were in 
good agreement with the benchmark result (example plotted in black). 

The examples shown highlight the importance of benchmarking to identify 
critical analysis input parameters. In summary, good agreement between the 
load-displacement relationship obtained from the propagation analysis results 
and the benchmark results could be achieved by selecting the appropriate input 
parameters. However, selecting the appropriate VCCT input parameters such 
as release tolerance, global or contact stabilization and viscous regularization, 
was not straightforward and often required an iterative procedure. The default 
setting for global stabilization yielded unsatisfactory results although the 
analysis converged. The detailed results for ABAQUS® are discussed in detail 
in a related report (Krueger, 2008). The results of a similar study using the 
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same approach for the finite element codes Marc and MSC.Nastran were 
published recently (Orifici and Krueger, 2010). 

5: CREATING A BENCHMARK SOLUTION FOR CYCLIC 

LOADING 

For the cyclic loading of the specimen, guidance was taken from a draft 
standard designed to determine mode I fatigue delamination propagation. In the 
draft document, it is recommended to start the test at a maximum displacement, 
dmax, which causes the energy release rate at the front, G/ max , to reach initially 
about 80% of Gic. The maximum load, P \naxi and maximum displacement, 
d max /2, were calculated using the known quadratic relationship between energy 
release rate and applied load or displacement given in equation (2). For the 
current study, a critical energy release rate G/ c =0.17 kJ/nr was used and the 
critical values P cr it and (%,u (grey dashed lines) were obtained from the 
benchmark for static delamination propagation shown in the load-displacement 
plot in Figure 6. The calculated maximum load, P max , and calculated maximum 
displacement, d max /2, are shown in Figure 6 (dashed red line) in relationship to 
the static benchmark case (solid grey circles and dashed grey line) mentioned 
above. During constant amplitude cyclic loading of a DCB specimen under 
displacement control, the applied maximum displacement, 6 max / 2=0.67 mm, is 
kept constant while the load drops as the delamination length increases (solid 
red circles and solid red line). The energy release rate corresponding to an 
applied maximum displacement 6 max /2=0.61 mm was calculated for different 
delamination lengths a using equation (2). The energy release rate decreases 
with increasing delamination length, a, as shown in Figure 7 (solid red circles 
and solid red line). Delamination growth was assumed to stop once the 
calculated energy release rate drops below the cutoff value, G t i„ (green solid 
horizontal line). The static benchmark case (solid grey circles and dashed grey 
line in Figure 7), where the delamination propagates at constant G/ c (solid blue 
line in Figure 7) was included for comparison. 

The number of cycles to delamination onset, No, can be obtained from the 
delamination onset curve which is a power law fit of experimental data 
obtained from a DCB test using the respective standard for delamination 
growth onset. Details are discussed in a related report (Krueger, 2010). Solving 
for Nd yields 


G = m 0 


■N 


D 



N d = c V G l \ 



( 3 ) 


At the beginning of the test, the specimen is loaded initially so that the energy 
release rate at the front, Gi max , reaches about 80% of G/ c corresponding to 
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6/^=0.1362 kJ/m 2 . From the delamination onset curve, the number of cycles 
to delamination onset is detennined as, No= \ 50. 

The number of cycles during stable delamination growth, Ng, can be obtained 
from the fatigue delamination propagation relationship (Paris Law). The 
delamination growth rate (solid purple line) can be expressed as a power law 
function 


da 

dN 


G n 


or 


A a 
A N 


G n 


( 4 ) 


where da/dN is the increase in delamination length per cycle and G max is the 
maximum energy release rate at the front at peak loading. The factor c and 
exponent n are obtained by fitting the curve to the experimental data obtained 
from DCB tests. A cutoff value, G>/„ was chosen below which delamination 
growth was assumed to stop. Details are discussed in a related report (Krueger, 
2010). Note, that this benchmarking exercise ignores branching, or fiber 
bridging. Hence, the Paris Law was not nonnalized with the static R-curve as 
recently suggested (Krueger, 2010). 

For the current study, increments of Aa= 0.1 mm were chosen. Starting at the 
initial delamination length ao= 30.5 mm, the energy release rates G umax were 
obtained for each increment, i, from the curve fit (solid red circles and solid red 
line) plotted in Figure 7. These energy release rate values were then used to 
obtain the increase in delamination length per cycle or growth rate Aa/AN from 
the Paris Law. The number of cycles during stable delamination growth, No, 
was calculated by summing the increments AN, 

k k ^ k 

N g = 2 A/V ' = ‘ Aa =* a = a 0 + y 2Aa = a 0 + k-Aa (5) 

1=1 ;=1 c ;-i 

where k is the number of increments. The resulting delamination length, a, was 
calculated by adding the incremental lengths Aa to the initial length ao 

For the combined case of delamination onset and growth, the total life, Nt, may 
be expressed as N T = N D + N c where, No, is the number of cycles to 
delamination onset and Ng, is the number of cycles during delamination 
growth. For this combined case, the delamination length, a, is plotted in 
Figure 8 for an increasing number of load cycles Nt. For the first No cycles, the 
delamination length remains constant (horizontal red line), followed by a 
growth section where - over Ng cycles - the delamination length increases 
following the Paris Law (crosses and solid red line). Once a delamination 
length is reached where the energy release rate drops below the assumed cutoff 
value, G t h, (as shown in Figure 7), the delamination growth no longer follows 
the Paris Law (dashed grey line) and stops (horizontal solid red line). A 
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delamination length prediction analysis that accounts for delamination fatigue 
onset as well as stable growth should yield results that closely resemble the 
plot in Figure 8. The curve fit (solid red line) can therefore be used as a 
benchmark. 

6: RESULTS FROM AUTOMATED GROWTH ANALYSES 

A set of example results are shown in Figure 9 where the delamination 
length, a, is plotted versus the number of cycles, N, for different input 
parameters and models. For all results shown, the analysis stopped when a 
10,000,000 cycles limit - used as input to tenninate the analysis - was reached. 

First, the effect of the initial time increment used in the analysis settings was 
studied as shown in Figure 9. The initial time increment was varied between 
4=0.01 (one tenth of a single loading cycle, 4=0.1 s, open blue circles and solid 
blue line) and 4=0.000 1 (one thousandth of a single loading cycle, open green 
diamonds and solid green line). For larger initial time increments, the onset of 
delamination shifted towards a lower number of cycles. Reducing the initial 
time increments, however, significantly increased the computation time. Based 
on the results, it was therefore decided to use an initial time increment of 
4=0.001 (open red squares and solid red line) for the remainder of the study to 
save computation time. This step is justified by the fact that the results obtained 
for 4=0.001 were almost identical to the values obtained from the analysis 
where a smaller initial time increment was used (4=0.0001). Second, the 
release tolerance was varied. In the current study, varying the input between 
0.2 (the default value) and 0.01 did not have any effect on the computed onset 
and growth behavior during the analysis. It is assumed that the release 
tolerance only affects static delamination propagation and not growth under 
cyclic loading studied here. Stabilization (to overcome the convergence 
problems) which was required for the automated delamination propagation 
under static loading, was not required for automated growth analysis under 
cyclic loading. However input parameters to define the cyclic loading were 
varied and showed different affects on the computed results. Also, solution 
controls had to be manipulated to achieve reasonable computation times. The 
results from these studies are discussed in detail in a related report (Krueger, 
2010). 

7: CONCLUSIONS 

The development of two benchmark examples for static delamination 
propagation and cyclic delamination growth prediction were presented and 
demonstrated for the commercial finite element code ABAQUS® Standard. The 
example was based on a finite element model of a Double Cantilever Beam 
(DCB) specimen, which is independent of the analysis software used and 
allows the assessment of the delamination propagation and growth prediction 
capabilities in commercial finite element codes. The development of the 



DEVELOPMENT OF BENCHMARK EXAMPLES FOR 
DELAMINATION ONSET AND FATIGUE GROWTH PREDICTION 


benchmark examples was presented step by step. Starting from an initially 
straight front, the delamination was then allowed to propagate or grow under 
cyclic loading using automated analyses. 

With respect to benchmarking, the results showed the following: 

• Benchmark examples for the assessment of automated delamination 
propagation and growth capabilities can easily be developed from a set 
of static analyses of models with different initial delamination lengths. 

• The development of benchmark examples is independent of the software 
used and independent of experimental results. 

• The benchmarking process can be used to identify the issues associated 
with the input of a particular finite element code or implementation 

• Benchmarking helps identify relevant input parameters so that they may 
be used with confidence when modelling more complex configurations. 

In general, good agreement between the results obtained from the propagation 
and growth analysis and the benchmark results could be achieved by selecting 
the appropriate input parameters. Overall, the results are promising. In a real 
case scenario, however, where the results are unknown, obtaining the right 
solution will remain challenging. Further studies are required which should 
include the assessment of the propagation capabilities in more complex mixed- 
mode specimens and on a structural level. 
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dimensions 



B 25.0 mm 

2/7 3.0 mm 

2L 150.0 mm 

a 0 30.5 mm 

layup: [0] 24 

static loading 

5/2 2.0 mm 

fatigue loading 

5 max /2 0.67 mm 
6 min /2 0.067 mm 
R 0.1 

f 10.0 Hz 


Figure 1: Double Cantilever Beam (DCB) Specimen. 



► x 

Figure 2: Deformed model of a DCB specimen (plane strain-CPE4 and plane stress-CPS4) 



Figure 3: Full three-dimensional (solid C3D8I) finite element model of a DCB specimen. 
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Figure 4: Calculated critical behavior and resulting benchmark case 
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Figure 5. Computed load-displacement behavior obtained from automated analysis 
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Figure 6. Critical load-displacement behavior for a DCB specimen. 
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Figure 7. Energy release rate - delamination length behavior for DCB specimen. 
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Figure 8. Delamination onset and growth behavior for DCB specimen. 
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Figure 9. Computed delamination onset and growth obtained for different initial time increments. 


